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In this paper we propose a systematic approach to construct mathematical models describing popu- 
lations of cancer-cells at different stages of disease development. The methodology we propose is 
based on stochastic Concurrent Constraint Programming, a flexible stochastic modelling language. 
The methodology is tested on (and partially motivated by) the study of prostate cancer. In particular, 
we prove how our method is suitable to systematically reconstruct different mathematical models 
of prostate cancer growth -together with interactions with different kinds of hormone therapy- at 
different levels of refinement. 
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1 Introduction 

The (long term) goal of this work is the development of a general framework for building programmable 
models to be used in the study of cancer. The models we have in mind are, in fact, models of populations 
of cells obeying specific growth and death laws and incorporating different stages of tumour evolution. 
At the architectural level, we envisage three key features: 

• programmability , as the possibility to enter a high-level description of the full model in terms of 
networks of interacting agents; 

• hybrid-ness, as the possibility to carry out a complete description of the semantics of our networks 
of agents by a dynamical system with a controlled combination of both discrete and continuous 
evolution; 

• stochasticity , as the possibility of specifying the interaction mechanisms obeying to given stochas- 
tic laws. 

The above mentioned networks naturally implement the hybrid nature of the different phases through 
which the development of tumoral populations undergo. Moreover, they allow us to study both the 
logical and quantitative features of a model, at a higher level of description than classical models based 
on Ordinary Differential Equations (ODE). As an example of such analyses, one can try to understand to 
what extent the observed noise depends on the structure of the internal/external interactions defining the 
model, as opposed to parameter variation. Finally, they allow to tune-up the level of discreteness to be 
used in the modelling activity. 
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Our concrete and motivating example is the study of prostate cancer, with special emphasis on the 
predictiveness ability of models including external interactions (in the form of medicament dispensa- 
tion) p0)[28| . For this reason we begin illustrating the various models of prostate cancer-cell evolution 
under different policies of chemical castration (see Section[2]and (28J). 

Our modelling approach to describe prostate cancer is based on stochastic Concurrent Constraint 
Programming (sCCP, [5]) a modelling language belonging to the class of stochastic Process Algebras 
(SPA) applied to Systems Biology fT2j . SPA have been applied in modelling many aspects of biological 
systems, including tumour growth (22j|23j. 

The objectives of this work are twofold: First we want to describe prostate cancer models and drug 
dispensation policies in sCCP. sCCP models of prostate tumour are presented in Section [4] Secondly, 
we want to clarify if noise observed in experimental data could be explained as a structural feature of the 
model. A first experimental study in this sense is presented in Section|5] 

sCCP is particularly suitable for the proposed task because its semantics is naturally stochastic, it 
is easily programmable and extensible, and it also has a general semantics in terms of stochastic hybrid 
automata [7] (see Sections 2.2 2.3 and 2.4 1. However, in order to properly use sCCP for cancer mod- 
elling, we had to enlarge the set of primitives that can be used to describe agent behaviour. In particular, 
we had to introduce a mechanism to describe actions triggered by conditions on model time and to give 
the agents the ability of changing their environment (i.e. system variables) according to random laws. 
These extensions turned out to be very simple to introduce in the hybrid semantics framework of sCCP 
(cf. Section [3). 



2 Background 

Below we briefly discuss the two main objects of our work in this paper: prostate cancer modelling 
techniques and sCCP, the modelling language that we will use to build our agent network in the rest of 
the paper. 



2.1 Prostate cancer modelling 

Prostate is a gland of the male reproductive system responsible for the production of seminal fluid. 
Prostate tumour is a very common (age related) disease, consisting mainly in the development of a mass 
of tumoral cells whose growth is correlated with (dependent from) the presence of androgen hormones 
(e.g. testosterone) [ fl0| . Most effective therapies consist in androgen deprivation (castration) by sur- 
gical or chemical means [ [TO] ]. Under such therapies, one observes an initial fast decrease of tumoral 
masses which, however, after a variable interval of time, undergo a relapse phase. This is caused by the 
emergence of a line of androgen independent tumoral cells, resistant to androgen deprivation |2|[T0). 

The above biological behaviours have been modelled using phenomenological approaches p9)20|28| 
based on ODEs. The variables of the model record the number of androgen dependent/independent cells 
and are equipped with growth and death laws expressing their time evolution. The spatial structure of the 
tumour is ignored, only the cell number is recorded. The observable of the model is the serum Prostate 
Specific Antigen p6| (PSA) — a bio-marker whose concentration in serum is strictly related to the size 
of the tumour mass, that is the number tumoral cells. 

More specifically, a model can be defined by the system of differential equations of Figure [T] which 



corresponds to the initial model used in [20] and [28 1. In the equations, we have three variables: x 



describes the population of Androgen Dependent (AD) cancer cells, y describes the population of An- 
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| = G x (x,z)-D x (x,z)-M xy ,(x,z) 
| = G y {y,z)-D y {y,z)+M xy {x,z) 
f = A-D z (z) 

P z = D z (z) = l 



G x {x,z) = a x (k 1 + (l-k 1 )^jAx 

D x (x,z)=fc(k 3 + (\-k 3 )^)x 
G y (y,z) = a y (l-dj 5 )y 
D y (y,z)=P y y 
M X} ,(x,z)=m l (l-f) 



Fi gure 1 ! ODE model of prostate tumor growth under Continuous Androgen Suppression, taken from |28| . P stands for 
Production, G for Growth, D for Death, and M for Mutation. In order to control P z the parameter a z is set to to represent the 
effect of androgen deprivation therapy. In absence of castration, such a rate would be P z = y. Parameters are set according 
to (28): a x = 0.0204, Ofy = 0.0242, /3 V = 0.0076, fi y = 0.0168, k x =0,k 2 = 2, k 3 = 8, Ar 4 = 0.5, mi = 5 ■ 10~ 5 , z = 20, d = 1, 
1 = 62.5. We agree that x controls the speed of z dynamics, zo is the stationary value of androgen hormone in absence of 
chemical castration, and rfe [0, 1] controls the effect of androgen hormone in the growth rate of AI cells (0 means that growth 
is independent from the hormone, 1 means that the growth is inhibited by the hormone). Notice how both the growth and death 
terms for AD-cells x depend non-linearly on the amount of hormone z through a sigmoid-shaped function 



drogen Independent (AI) cancer cells, and z is the concentration of androgen hormone. The concentration 
of PSA is computed simply as x + y. This model is able to capture the relapse phase of the tumour due 
to androgen independent cells growth 1 19 20] . 

A (clever) clinical approach tackling the problem of relapse after a long interval of chemical cas- 
tration, is the so-called Intermittent Androgen Suppression (IAS) policy, see (2j. The effectiveness of 
the IAS therapy is based on the fact that the underlying mechanism of the disease seems to involve a 
competition between androgen dependent and androgen independent cells. This implies that a complete 
reduction of the number of AD cells removes any obstacle to AI cells growth. Hence, a more effective 
strategy consists in maintaining a certain number of AD cells to inhibit AI cells growth. The overall 
effect of an intermittent androgen deprivation strategy is the delay — possibly for a very long time — of 
the development of an androgen independent tumour. 

Given the quantities involved and the presence of an external input — the medicament dispensation 
strategy — to keep into account as a time-controlled phase change, the mathematical modelling machinery 
naturally evolved into a hybrid model. To be precise, the system can switch between states representing 
normal and androgen deprivation modes [19 28 1. These modes can be conveniently represented by a 
boolean variable u, where u = 1 describes the drug dispensation phase. The equations are obtained from 
the set of Figure[T]by replacing P z in the equation for the androgen hormone z by the w-dependent function 
P z (u) = Note that P z (\) = 0, as in the previous model. A further evolution of the model consisted 

in the introduction of noise in the equation, thereby moving to a set of stochastic differential equations. 
This move was performed as an attempt to capture small fluctuations observed in the PSA data. 

Our first motivation was to attempt to clarify whether the noise introduced in the model of Aihara et 
al. p8[ is external or internal. In the former case, the explanation would call into play additional (noise) 
sources not related with growth and death laws of AD/AI cells. In the latter case, the noise could be 
explained in terms of fluctuations of such transitions, when considered as stochastic. Our approach is 
to begin designing a discrete and stochastic version of the above model, in terms of Continuous Time 
Markov Chains (CTMC) |l7|[25j. The construction will be given in Section [4] and the experimental 
results will be given in Section [5] Below we briefly describe our agent's language, that will act both as 
an intermediate layer in the translation to CTMC and as a computational counterpart of the model based 
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on (stochastic) differential equations. 



2.2 Stochastic Concurrent Constraint Programming 

We briefly introduce now (a simplified version of) stochastic Concurrent Constraint Programming (sCCP), 
sketching the basic notions needed in the rest of the paper. More details on the language can be found 
in Q. sCCP has two basic ingredients: agents and constraints. Agents are the main actors, interacting 
by asynchronously exchanging information, in form of constraints, through the constraint store. sCCP 
has been mainly applied as a modelling language for biological systems BJ, using the constraint store 
to describe the state of the system, e.g. numerosity of molecular species. These quantities are described 
by a set of variables that can change value during computation, called stream variables [5]. At least 
for modelling simple biological scenarios, one needs very simple constraints, basically comparing and 
assigning new values to stream variables. 

Definition 2.1. An sCCP program is a tuple 21 = (A,D ,X,init(X)), where 

1. The initial network of agents A and the set of definitions D are given by the following grammar: 

A =M | A || A M = n.C\M+M n = [G(X) -> R(X,X')] X (x) 



D = | DUD | {C d =M} 



2. X is the set of stream variables of the store (with global scope), usually taking integer values; 

3. initiX) is a predicate on X of the form X = xo, assigning an initial value to store variables. 

In the previous definition, basic actions are guarded updates of (some of the) variables: G(X) is a 
quantifier-free first order formula whose atoms are inequality predicates on variables X and R(X,X') is a 
predicate on X, X' of the form X' = X + k (X' denotes variables of X after the update^] for some vector 
k £ W. Each such action has a stochastic duration, specified by associating an exponentially distributed 
random variable to actions, whose rate depends on the state of the system through a function A, with 



values A(X) £ M + . The semantics of sCCP is given by a Continuous Time Markov Chain [25], deduced 



from the labelled transition system associated with an agent netw ork by the operational semantics (5J. 



def 

All sCCP agents C = ' M £ D, according to Definition 



2.1 



are sequential, i.e. they do not contain 



any occurrence of the parallel operator, whose usage is restricted at the upper level of the network. sCCP 
sequential agents can be seen as automata synchronizing on store variables and they can be conveniently 
represented as labelled graphs, called Reduced Transition Systems (RTS) (see [6]). More precisely, 
RTS(C) = (S(C),E(C)) is a multi-graph with vertices 5(C) corresponding to the different states of an 
agent and with edges e £ E(C) corresponding to actions performable by the agent, labelled by the cor- 
responding rate (k e ), guard (G e ), and reset (R e ). In the following, we also need the notion of extended 
sCCP program 2l + , in which we introduce a new set of variables P = {Pc \ C £ D} taking integer values 
and recording how many copies of agent C £ D are in parallel in the current sCCP-network. Rates, 
guards, and resets are modified to consistently treat P-variables. With this trick, we can assume that in 
2l + there is never more than one copy of the same agent running in parallel at the same time. Further 
details on these notions can be found in |6j. 

Examples of sCCP-programs will be given below and throughout the rest of the paper. We begin by 
observing that playing with the meaning assigned to variables and transitions, models at different levels 



'The constraints that can be used to update the constraint store are rather limited, as they simply add a constant to some 
stream variables. This restriction, however, allows to interpret sCCP-actions as continuous fluxes, a required condition to define 



the hybrid semantics, see also Section 2.4 
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of detail can be programmed in a multi-level environment. At the lowest level, even a single cell can be 
represented by an agent, thereby suggesting a programming style such as: 

cell_H(M) :- [M <Ho^M' =M+ l] A(i ( M) .cell_H(M) + [M = Hq -> M 1 = M] Xm .cell_T(M) 
cell_T(M) :- [M < -> M' = M+ % (M) .cell_T(M) + [M = -> M' = M] v 

In the previous piece of code, we are describing the mutation from a cell in a healthy state (cell_H) to a 
cell in a tumoral state (cell_T), as a consequence of accumulation of (internal) genomic mutations. The 
number of mutations of a cell is represented by a (local) variable M, which is increased by one when a 
mutation event occurs. A cell will remain in a healthy state as long as the total number of mutations is 
bounded by /lo (first action). We assumed that mutations can happen at a rate dependent on the internal 
state of the cell (number of mutations M). When such a critical level is reached, the cell after some time 
will become tumoral (second action). A tumoral cell, instead, will survive while the number of genomic 
mutations remains below pt\ (third action). When such a threshold is reached, the cell will die (fourth 
action). A dead cell is represented by the null agent 0, which is removed from the system. 

The full program would run in parallel a number of agents equal to the number of cells modelled, 
each with its specific mutation-counting variable. We remark that this very simple model, not including 
events like cell duplication and cell apoptosis, has been introduced just for illustrative purposes. 

The above (low) level of modelling will not be used in the following, even though along the lines of 
the above example the competition mechanism between AD and AI cells could be explicitly programmed. 
Instead, we will use variables counting the number of AD and AI cells, and assign agents to interactions, 
like cell population growth/death. This corresponds to "programming the dynamics" at a higher level of 
abstraction. We will discuss this usage in the next section. 

Even if the standard semantics of an sCCP-program is given in terms of CTMC, sCCP-programs 
have also different semantics, based on ODE [6] and hybrid automata [7], providing an additional degree 
of flexibility to the language. In particular, the hybrid semantics is parametric with respect to the degree 
of continuity introduced in the model, and it subsumes as special cases both the CTMC and the ODE 
semantics. We will give now more details on this semantics, exploiting it in Section[3]to smoothly extend 
sCCP with additional primitives that we will need for modelling tumour-cells]^] 

2.3 Transition-Driven Stochastic Hybrid Automata 

Transition-Driven Stochastic Hybrid Automata (TDSHA) have been introduced in (7j|9j as a convenient 
formalism to associate a stochastic hybrid automaton to an sCCP program. 

The emphasis of TDSHA is on transitions which, as always in hybrid automata, can be either discrete 
(corresponding to jumps) or continuous (representing flows acting on system's variables). More specifi- 
cally, there are two kinds of discrete transitions: 

1 . instantaneous or forced: taking place as soon as their guard becomes true, and 

2. stochastic: occurring after an exponentially distributed delay. 

Definition 2.2. A Transition-Driven Stochastic Hybrid Automaton (TDSHA) is a tuple 
T= (e,X,^f,y,mif), where: 

• Q is a finite set of control modes and X = {X\ , . . . ,X n } is a set of real valued system's variable^ 

2 A software tool to model and analyse sCCP programs is under development. A preliminary version is available from the 
authors upon request. 

3 Notation: the time derivative of Xj is denoted by Xj, while the value of Xj after a change of mode is indicated by Xj 
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• ^ is the set of continuous transitions or flows, containing triples z = (q,s,f), where: q G Q is a 
mode, s is a real vector of size |X|, and / : W — > M is a (Lipschitz) function. They are indicated 
by q T , s T , and f x , respectively. 

• @ is the set of discrete or instantaneous transitions, whose elements are tuples r\ = (q\,q2,G,R,w), 
where: q\ is the exit-mode, q2 is the enter-mode, and w : K n — > M + is a weight function used to 
resolve non-determinism between two or more active transitions. Moreover, G is a quantifier-free 
first-order formula with free variables among X and R, which is a conjunction of atoms of the form 
X' = r(X,/i), where r : W +m — > M, is the reset function of X, depending on system's variables X 
as well as on a vector of parameters 11. In particular, parameters in 11 can be either constants or 
random values. The elements of a tuple r\ are indicated by q* , , Wf, , , and , respectively. 

• 5? is the set of stochastic transitions, whose elements are tuples r\ = (q\,q2,G,R,X), where q\, 
q2, G, and R axe as for transitions in while X : K" — > M + is a function giving the state-dependent 
rate of the transition. Such function is referred to by Xq. 

• init is a pair (go,Xo) G Q x W, identifying the initial state of the system. 

The dynamics of TDSHA can be formally defined (7J by associating with them a Piecewise Deter- 
ministic Markov Processes [14]. Intuitively, the TDSHA dynamics is given by periods of continuous 
evolution, interleaved by discrete jumps, as customary with hybrid systems. 

Within each mode q G Q, the TDSHA evolves following the solution of a set of Ordinary Differential 
Equations (ODE), constructed combining the different continuous transitions active in q. Each such 
contributes to the ODE with a flow given by f T (X), whose effect on each variable is described by 
s T . Hence, in mode q G Q, ODE are given by X = £ Te <^ \q\= q s T f t(X). 

Two kinds of discrete jumps are possible: stochastic transitions rj G 5? are fired according to their rate, 
while instantaneous transitions r\ g <2i are fired as soon as their guard becomes true. In both cases, the 
state of the system is reset according to the policy specified by R, 7 . Choice among several active stochas- 
tic or instantaneous transitions is resolved probabilistically according to their rate or priority. 
We stress that TDSHA have another sources of stochasticity in addition to priorities and rates: resets. 
Resets function, in fact, can assign random values to variables, through their dependency on random 
variables (via the parameter's tuple This feature of TDSHA will be exploited in Section[3] 
Furthermore, a notion of asynchronous TDSHA product T = Ti ® %2 can be defined in a simple way, 
see (9j. Essentially, the discrete state space of the product automaton is Q\ x Q2, while transitions from 
state (gi,^) are all those issuing from q\ or q2. 

2.4 Hybrid Semantics of sCCP 

In this section we briefly recall the definition of the hybrid semantics for sCCP |7J, which is given by as- 
sociating a TDSHA to an sCCP program. The mapping is compositional: first single sCCP components 
are converted to TDSHAs, then one TDSHA is obtained as their product. 

A TDSHA is, ultimately, an automaton and the discrete skeleton of the TDSHA associated to an 
sCCP component C is directly derived from the RTS -which is a labelled graph- of C. The level 
of discreteness/continuity will be a parameter, which is fixed by partitioning RTS -edges E(C) into 
discrete edges E^C) and continuous edges E C (C). The former ones will generate stochastic tran- 
sitions and the latter ones continuous transitions. The edge partition is encoded by a boolean vec- 
tor K of size E(C), where K[e] = 1 means that e G E C (C). The recipe for constructing the TDSHA 
^c, k = ( Qc > Y , j j 1 initc ) of component C of 2t + is sketched below. 
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Discrete Modes. Continuous edges in E C (C) can short-circuit different RTS-states of S(C). Hence, we 
collapse together RTS-states connected by an (undirected) path of continuous edges. Let [Q] c be the 
set — equivalence class — of states containing C, E S(C). Then Qc = {[Q] c | C, E 5(C)}. 
Variables. The TDSHA uses exactly the variables Y of the extension 2l + : system variables X and state 
variables P. 

Continuous transitions. Each edge in e E E C (C) becomes a transition in (q,s,f) E ^c> where: q = [Ci\ c , 
with C, the source state of e, s is the update vector defining R e (i.e. R e equals Y' = Y + s), and f(Y) = 
A e (Y). 

Stochastic transitions. Each edge e E E C /(C) produces a discrete stochastic transition (q\,q2,G,R,X) E 
J^c, where q\ is the equivalence class containing the source state C\ of e, q2 contains the target state C2 

of e, G = G e , X = X e , and R = R e Ixstatevarjres. The second term of R is needed to correctly deal 

def 

with classes of states of S(C): statevar_res = (Ace^i P'c = 0) ^^c = 1- 

Initial conditions, initc is constructed by considering the initial state of the sCCP component in 2l + . 
Instantaneous transitions. Instantaneous transitions are not needed at this level of the mapping. They 
have been used in [7 ] to describe dynamic partition policies between discrete and continuous transitions, 
and they will be used in the next section to give sCCP new functionalities. 



Once we have constructed the TDSHA for each parallel component of 21+ corresponding to the initial 
network A = C\ \\ ... \\ C n , we apply the (basically) standard product construction introduced at the end 
of|23]to obtain Ta,k,©...©k„ = Z Cl , Kl ® ■ • -®^c n ,K n - 



Lattice of TDSHA. The reader can easily see that the previous definition is parametric with respect 
to the degree of continuity introduced in the hybrid automata. We can arrange the different TDSHA 
obtained by different choices of K in a lattice, where the top element is the fully continuous and deter- 
ministic TDSHA, while the bottom element is the fully discrete and stochastic TDSHA. In particular, 
the fully continuous TDSHA has one single mode, no discrete transitions, and it corresponds to the set 
of ODE associated to an sCCP program by fluid-flow approximation (6j. On the other hand, the fully 
discrete TDSHA is the CTMC associated with an sCCP by the standard semantics. 

We stress the fact that sCCP is not a language designed to describe stochastic hybrid systems, but it 
is rather a programming language to model networks of (stochastically) interacting agents. The hybrid 
semantics has been defined as a computationally efficient approximation of the stochastic one. Con- 
sequently, the sCCP syntax has no way of specifying if an action has to be interpreted as discrete or 
continuous. This choice has to be performed at the RTS level. Practically, we are investigating the def- 
inition of general rules providing good automatic partitions (in terms of behaviour approximation and 
computational efficiency) and the external annotation of sCCP actions by naming them, similarly to Bio- 
PEPA fT3| . As a final remark, we point out that the asynchronous nature of sCCP allows the definition 
of the hybrid semantics in a relatively simple compositional way. 



3 Extending sCCP: Events and Random Updates 

In this section we will extend the sCCP language by equipping it with additional primitives that will be 
used, in our case study, to model prostate tumour dynamics. In particular, in order to model drug delivery 
policies, we need to describe SBML-like (T]|TT| events in sCCP. Furthermore, in the experiments we 
will carry out in Section |5j we will also need to reset variables to random values, drawn according to 
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given distributions. In the following, we will briefly introduce these extensions, both syntactically and 
semantically. More specifically, we will add new syntactic primitives to sCCP, framing their semantics 
within the skeleton of TDSHA, extending the hybrid semantics discussed in the previous section. 



Events. SMBL-like events describe instantaneous modifications of the system, triggered by conditions 
on the system state or on the simulation time. Events allow us to model very easily external influences 
on a system, such as, for instance, the hormone deprivation policies for prostate cancer. More precisely, 
following pT[ , an SBML-like event is a conditional expression of the form 

if (trigger) then (assignments) with (delay), 

where trigger is a condition involving system variables and the simulation time, assignment is an 
update of (some) system variables, and delay is an optional time delay between the activation of the 
trigger and the consequent update of the variables. 

We will first discuss how to describe in sCCP events with trigger depending system's state only 
and with no delay. In order to do this, we simply have to extend sCCP with instantaneous actions, i.e. 
actions taking no time to occur. Syntactically, this can be done by allowing rates to take value infinity, 
i.e. we have to allow basic actions of the form [G — > /?]«. Even if instantaneous transitions are quite 
standard in stochastic process algebras (3 15 T8j, we will take care to define their semantics within the 



TDSHA framework. In particular, we will map an action C = [G — > R]„.Ci +M to the instantaneous 
transition ([C] c , [C\] c ,G,R,w), where the priority w is the constant function 1 (i.e. w(X) = 1 for each X). 
Therefore, these sCCP actions are always kept discrete while mapping sCCP-programs to TDSHA. If we 
consider the TDSHA obtained by keeping all sCCP transitions discrete, then the sCCP semantics with 
instantaneous actions remains a CTMC, provided the sCCP program cannot execute an infinite number 
of instantaneous actions in the same time instant (a well-known fact, cf. [15]), because non-determinism 
between two or more active instantaneous transitions of TDSHA is resolved probabilistically by the use 
of priorities. 

Time controlled events, instead, are more difficult to express in stochastic process algebras, as they 



change the semantics of the language, from a CTMC to a semi-Markov process [24 1, i.e. a stochastic 
process in which transition times can be sampled according to general distributions. In particular, we 
are considering semi-Markov processes mixing exponential and deterministic transition times. Such 
processes can be simulated using priority queues, as customary done in discrete event simulation. The 
introduction of time controlled events in sCCP, however, is done by simply adding a special (reserved) 
variable, Time, to the language. The use of variable Time is restricted to infinite-rate actions, like [Time = 
10 — > R]cc, and it cannot be reset. 

From a semantic point of view, the value of Time corresponds to the simulation time of the model. 
This fact can be built into the TDSHA framework: we add Time to the TDSHAs variables and we 
render its dynamics by adding a new automaton (sometimes called time-monitor) in parallel with the 
TDSHA obtained from the sCCP program. This automaton will have one single mode q, one variable 
(Time), initial state (q,0), and one single continuous transition, (q,s,f) — modeling the flow of time 
with rate / = 1 and stoichiometric vector s = 1. As all actions containing Time have infinite rate, they 
will all become instantaneous transitions in the TDSHA. Therefore, an action with guard Time = 10 
will activate as soon as TDSHA variable Time, which corresponds to the model time, will reach value 
10. Of course, the guards of infinite-rate actions can combine in complex ways conditions on time and 
on system variables. Notice that events with a delay d between the activation of the trigger and the 
consequent assignment of system's variables can be easily modelled within sCCP as a sequence of two 
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instantaneous actions: [trigger — > T = Time + d]cc. [Time = T — > assignment}^. We can formally render 
the correctness of this construction in the following 

Theorem 3.1. The TDSHA obtained from an sCCP program with instantaneous transitions and time 
constraints, keeping all the non-instantaneous transitions stochastic and coupling it with a time-monitor, 
is stochastically equivalent to a semi-Markov process. 



Random Updates. The capability of updating system variables with random values can be useful in 
many applications, and it will be exploited in this paper to model uncertainty and intrinsic variability in 
(kinetic) parameters. For instance, we can model an action whose rate is a random variable rather than a 
number, abstracting in this way an (unknown) underlying mechanism giving rise to such randomness. 

Alternatively, random resets can be used to incorporate uncertainty about a parameter in the model, 
according to the Bayesian point of view |29|. In particular, we can select randomly the value of a 
parameter at the beginning of each simulation according to a given prior distribution (reflecting our 
knowledge on the parameter value). Then, the empirical distribution of the certain and the uncertain 
model are compared, assessing the sensitivity of the model with respect to the parameter^] 

To formally add random resets in sCCP, we allow arbitrary reset functions in infinite-rate sCCP 
actions [G — > /?]«,. Therefore, the reset R of such an actions can now be a conjunction of atoms of the 
form X' = /?(X,P,ju), where pt is a vector of parameters, which can be constants or random variables, 
like in the definition of TDSHA. The mapping from sCCP to TDSHA does not need to be modified: the 
instantaneous TDSHA transition associated to an infinite-rate action is the one described for events. We 
observe that there is no real obstacle in using general resets also in sCCP actions with finite rate, but in 
this case these actions cannot be approximated continuously in the hybrid semantics, but need to be kept 
discrete. 

We stress that the advantage of defining the semantics of these extensions in terms of TDSHA (in- 
stead of defining them at the operational semantics level, which is possible at least for random resets) 
is that they apply to all semantics of sCCP, including the CTMC and ODE based semantics, which are 
special instances of the TDSHA semantics. 

Equipped with these extensions to sCCP, we now focus on prostate tumour growth modelling. Our 
first goal is to build a discrete and stochastic (programmable) model, capable of rendering the prostate 



tumour-growth model(s) presented in Section 2.1 This approach will be discussed in the following 



section, while the effects of noise in this model will be discussed in Section [5] 



4 sCCP model of prostate cancer growth 

In this section we describe how to go from a mathematical model like the one of Section [2] to sCCP- 
programs, illustrating the technique directly on the prostate tumour case. 

The basic principles in our approach consist in identifying from the differential equations the sCCP 
variables and the sCCP interactions. Variables will roughly correspond to variables in the differential 
equation model. Interactions will be obtained by disassembling the right-end sides of the differential 
equations, identifying explicitly the different actions modifying the populations. Interactions will be 
described by agents of the network. 

For prostate tumour cancer, following | [28j , we use four variables X,Y,Z, and V standing for the 
(numbers of) AD cells, AI cells, androgen hormone molecules, and PSA molecules, respectively. 

4 An alternative way to incorporate lack of knowledge is to use imprecise probabilities, like in Imprecise Markov Chains (271. 
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growthAD :- [X > -> X' = X + l] Gx ( X Z ).growthAD deathAD :- [X > -¥ X' = X - i] Dx r x , Z ).deathAD 

growthAI :-[Y >0^Y' = y + l] Gf ,( FZ) .growthAI deathAI :- [Y > -> F' = y - l] Dl , (FZ) .deathAI 

produceANDH :- [?rae -> Z' = Z+ %.produceANDH degradeANDH :- [Z > -> Z' = Z- l] Dz(z) .degradeANDH 
producePSA :- [true -*V' = V + 1]/\,(X -producePSA degradePSA :- [V > -* V' = V - l^^.degradePSA 
mutateADtoAI :- [X > -» X' = X - 1 A Y' = F + . Z ).mutateADtoAI 

Gx(X,Z) = a v ^,+(l-^ 1 ) ZT | I¥ U Z) y (y,Z) = /3 y y Dz(Z) = f 

Jfc(X,Z)=&(* 3 + (l-* 3 ) z 4 I5 )x MxF(Z,Z)=m 1 (l-^) /V(X,y) = ^<|±^ 

G F (y,z) = a v (i-rf a ^)y p z = o d v (v) = v 

Figure 2: List of agents of the sCCP-program modelling the prostate tumour growth. The initial configuration of the network 
consists of all the agents above, running in parallel. Rates are obtained from the ODE of Section[2]by suitably scaling parameters 
according to conversion factors, to be denoted fl z , £ly, and Nq. £ly is the conversion factor between PSA concentration and 
PSA numerosity. It corresponds to the approximate number of molecules in a nanogram per millilitre of PSA. Analogously, 
Q.Z is the conversion factor for Z, i.e the number of molecules giving a concentration of a nano-mole per litre. A'o, instead, is 
the reference number of cells, defined as the number of cells that on average produce 1 3j? units of PSA. The production and 
degradation rates of PSA are defined so that the average stationary number of PSA molecules is corresponding to a 

concentration equal to x + y, which is the PSA value computed in the ODE model of Section|2] 

It will be convenient to classify interactions into two classes: cellular interactions (i.e. growth, death, 
mutation) and molecular interactions (i.e. production, degradation). For prostate tumour cancer we 
have five cellular interactions, corresponding to growth and death for AI and AD cells, in addition to an 
interaction modelling mutation of AD into AI cells. As far as molecular interactions are concerned, we 
consider four of them, that is production and degradation of androgen hormone and PSA. The differential 
equation model defines the PSA level as the sum of the number of AI and AD cells. We have chosen, 
instead, to treat PSA production and degradation explicitly, in order to give more internal flexibility to 
the model. 

In Figure [2] the reader can find the full network of agents corresponding the equations presented in 
Section [2] In order to illustrate how agents interact, a precise definition of rates must be given. Consider, 
for example, the agent corresponding to growth of AD cells: 

growthAD :- [X > -> X' = X + l] Gx(XZ ). growthAD 

The agent represents a loop in which the number of AD cells grows by one at each iteration. The growth 
takes place only when the guard X > is satisfied, and it has the effect of incrementing the value of X 
by one unit. The rate of the interaction is directly derived from the differential equation models: 

Gx(X,Z) = o x (k l ^-h)z^)x. 

The parameter Q.z appearing in the above rate, is not found in the differential equation model: it is a 
scaling parameter necessary to convert the concentration z of the hormones into the molecular count 
Z. Other conversion factors are required for the size of cell populations X and Y (No) and for PSA 
concentration (Hy), see caption of Figure[2] The significance of these conversion factor in the dynamics 
of the entire model will be discussed in the Section [5] 

4.1 Hybrid dynamics and drug dispensation 



The sCCP-program described above corresponds to the Continuous Androgen Suppression (CAS) pol- 
icy [10 20 1 . As said, more effective drug dispensation policies have been studied and in particular the 
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(b) ODE vs Stochastic, IAS policy 
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Figure 3: 1 3(a) and 3(b) i Comparison of simulated trajectories of PSA in the sCCP-programs and in the ODE models, 
both under CAS and IAS drug dispensation policies. Stochastic trajectories are generated by a Gillespie-like simulation algo- 
rithm (17| . The time unit is a day. Parameters of ODE are as in Figure[TJ Scaling parameters of the stochastic model are set to 
N{) = 10 s , Q.z = 10 5 , fly = 10 5 . PSA of the stochastic model is rescaled as v = ^- before being plotted, so to have the same 
scale in both modes. Initial values of scaled variables are x= 15, y = 0.1, Z= 12. In the CAS policy, chemical castration is 
always in force. The IAS policy adopted is the following: every 4 weeks (28 time units) the value of PSA is checked. Drug 
dispensation is interrupted if normalized PSA has dropped below 4, while it is resumed when PSA exceeds 10. The behaviour 
of the stochastic model is essentially indistinguishable from the behavior of the ODE one. This basically happens for scaling 
parameters greater than 10 3 , hence their precise values do not have a relevant impact in the dynamics, as long as they are large 
enough. ( |3(c)| > Comparison of single trajectories for CAS policy, and different values of Nq. %3(d)\ Single trajectory for IAS 
policy in the presence of an external random disturbance source. 



Intermittent Androgen Suppression (IAS) has been proven effective to control prostate tumour develop- 
ment [2 1. The mathematical rendering of such a policy calls naturally into play a hybrid model describing 
the on/off modes of drug dispensation. In our stochastic program this is done by introducing a variable U 
(corresponding to the variable u of the differential equation model) governing the switch between on/off 
androgen deprivation policy. The syntactic feature of U is that it turns out to appear only in guards and 
not in rates: it is purely a control variable. More specifically, it appears in the new agent corresponding 
to androgen hormone production: 
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produceANDHc :- [U = -» Z' = Z + 1] - Q n z .produceANDHc + [U = 1 -)• Z' = Z+ l] .produceANDHc 

The above agent replaces the corresponding agent (produceANDH) in Figure [2] 

Typical drug dispensation policies for androgen deprivation, control PSA concentration at fixed time 
intervals (usually every 4 weeks) and determine whether dispensation should be resumed/suspended by 
checking whether specific threshold values are reached. In order to describe such a policy, we can use 
sCCP agents with infinite rate and guards on system time, as described in Section [3j More precisely, we 
just need a variable W, which will describe in which day the level of PSA will be controlled, initially set 
to 8 t = 28, and add in parallel to the sCCP program the agent checkPSA.on, defined by: 

W + 8, A U' = 0]oc.checkPSA_of f 
W + c5,]oo.clieckPSA_on 
= W + 8 t AU' = l]oo.checkPSA_on 
= W + (5 r ]oo.checkPSA_of f 

Notice that, according to the semantics of Section [3j this is a proper hybrid model, as the explicit 
management of simulation time is a continuous ingredient. However, our semantics machinery in terms 
of TDSHA allows a much higher degree of flexibility, for instance treating (some of) the system variables 
as continuous to speed up simulation times. 



checkPSA.on :- [Time = W A V < v on -> W = 

+ [Time = W A V > v on -> W' = 

checkPSA_of f :- [Time = W A V > v off -> W' = 

+ [Time = W A V < v // -> W' = 



5 Experimental results 

In this section we briefly present an experimental analysis of the models of the previous section. More 
details on the analysis can be found in the supplementary material, available online EJ. The parameters 
are fixed to the same values of the ODE model (see Figure[T]). The three scaling factors, A^o, O-z and Cly, 
have been fixed to values that we deemed meaningful and that were checked against real data, see the 
caption of Figure 



' 1 



The first analysis that we consider is the comparison of the evolution of the agent-based stochastic 



model with the temporal evolution of ODE. In Figures 3(a) and 3(b) we compare the dynamics of PSA 
for the CAS and IAS drug dispensation policies, respectively. We show only the value of PSA, as it is 
the only observable quantity of the model. As we can see, the two models behave essentially in the same 
way. Even if a similar behaviour had to be expected, the fact that the stochastic system basically shows no 
noise at all may be seen as quite surprising. Actually, this phenomenon can be easily explained observing 
that the variables we are considering are all taking large values, i.e. they correspond to large populations, 
on the order of millions of cells or millions of molecules (in Figure [3] a PSA level of 10 corresponds 
to 1 million molecules). In these circumstances, the relative magnitude of fluctuations, which is of the 
order of is too small to produce significant effects 1 161. This essentially means that the variability in 
behaviour between single cells is lost when we consider large populations: the differences cancel out and 
the observed behaviour essentially coincides with the average one (cf. also Section 1 of Supplementary 
Material gj). 

One of the advantages of having a discrete and stochastic model is that we can also study the be- 
haviour of the model when the population of cells is small. One way to perform this experiment in this 



5 Changing the value of the scaling factors can be seen as assuming a different granularity in counting. For instance, if we 
have N tumour cells, and we change the reference parameter Nq from 1 to 10, it means that we count how many groups of 10 
cells are present in the system. 
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setting is to decrement the parameter No, with the overall effect of reducing the total number of cancer 



cells. In Figure 3(d) we show the evolution of the number of tumour cells for different values of No. As 
expected, as No decreases the evolution becomes more noisy. Interestingly, for very small populations 
the model changes its behaviour: there is a considerable probability that AI cells can extinguish before 
AD cells are eliminated by chemical castration, so that the tumour tends to be completely eliminated and 
no relapse can be observed. These low numbers may describe a situation during the initial stages of the 
tumour. Interestingly, the IAS therapy is less effective for low populations, as it reduces the chance of 
extinction of AD cells, so that the probability of tumour extinction is considearably lower. Details can 
be found in Section 2 of the Supplementary Material fl4j. 

The sCCP-program model of tumor growth of Section [4] essentially does not show any noise in case 
of large cellular populations. Therefore, we considered possible modifications of the model to introduce 
some form of internal variability. This approach can be justified in order to check whether the variability 
in PSA concentration that is observed in real measurements, can be explained by the simple structure 
of this model. First, we modified the model adding a random source of variation in the PSA production 
rate. In Section 3 of the Supplementary Material [4] we present an experiment in which the production 
rate of PSA is no longer a constant with respect to the total number of tumor cells, but it is variable. 
In particular, we assume that production rate of PSA is a random number, uniformly distributed in the 
interval [0,2 ^ v ^ +r ^ ]. This is accomplished in the sCCP-program by introducing a new variable K (the 
new production rate) and replacing the agent producePSA with the following one: 

producePSA' :- [true ->• K' = UniffO^U.ltrue ->■ V = V + l]^a v fv v , .producePSA'. 

Interestingly, even a large randomness in the PSA production rate does not result in a significant noisy 
behaviour: also in this case the fluctuations are averaged out. 

We also considered a different modification of the model, in which we try to see whether noise may 
emerge as a consequence of the interaction with some external mechanism, possibly dependent on some 
global (physiological) condition. Specifically, we added a variable H (for "hidden") that can assume just 
two possible values, namely and 1, representing the presence/absence of an unspecified (physiological) 
condition. We assume that H switches from to 1 and vice versa twice a day on average, and that the 
production rate of PSA is subject to a threefold increase when H = 1. Hence, we added to the model the 
following agent: 

hidden :- [H = -»• H' = 1] 2 + [H = 1 ->■ H' = 0] 2 , 
and we replaced the PSA degradation rate by (1 +2H). Trajectories for this model are shown in Figure 



3(d) In this case, we observe noise also in the case of large populations, cf. Section 4 of the Supplemen- 
tary Material (4}. 

We stress that these last two models do not necessarily have any biologically significant interpreta- 
tion. They are just illustrative examples of possible sources of noise in systems with large populations. 

The analysis we carried out is different from the one of [28]. In this paper, noise is introduced in the 
model as an additional disturbance term in the ODEs (hence using Stochastic Differential Equations). In 
this case, however, the model does not provide any possible explanation of noise in terms of mechanisms 
intrinsic to the system under study. On the contrary, our study is aimed at better clarifying if the observed 
noise can be explained in terms of specific mechanisms, i.e. in terms of intrinsic properties of a given 
model. Notice that some noise will inevitably be introduced by factors external to the system, for instance 
by measurement errors. 

What we understood from the analysis presented here is that the structure of phenomenological 
tumour cell growth models, like the one considered in this paper, may not be sufficiently rich to contain 
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internal mechanisms for noise generation. If one is interested in these issues, then more complex models, 
taking into account more detailed biological mechanisms, should be considered. These models can also 
be easily described in our programming framework. We stress that this kind of analysis is more easily 
carried out in a discrete and stochastic setting. 

6 Conclusion 

The technique presented in this paper consists in showing how to step from a differential equation model 
to a "program" model, taking the form of a network of interacting agents. The specific programming 
language we used, allows us also to introduce a stochastic element (internal to the model and) rendered 
as the speed at which any specific interaction takes place. In particular, we considered a few extensions 
of sCCP-programs to describe time-driven events and random updates. General sCCP-programs, in ad- 
dition, can easily model cell duplication events in the low level model of Section [2] by an unrestricted 
usage of parallel composition and local variable declaration. In general, communication between agents 
representing single cells can be easily encoded in an asynchronous setting like the one of sCCP using 
dedicated variables, playing the role of communication channels, or modelling protein-mediated interac- 
tion. Exploiting the programmability of the shared memory (constraint store), one can easily introduce 
spatial information [ 8 ] or more complex cell interaction rules. Hence, sCCP-programs allow us to model 
explicitly geometrically qualified interactions or complex competitive dynamics regulating cells growths 
and deaths. However, a satisfactory definition of the hybrid semantics for this larger class of sCCP 
programs is still an open issue. 

A programming environment like sCCP allows the construction of a "wizard" for fast prototyping of 
(cancer) cell population dynamics. Among other things, this approach should allow us to easily address 
such basic questions as the effect and nature of noise, parameter dependencies, logical structure of the 
interactions, etc. More advanced analysis techniques, like statistical model checking (2TJ , will further 
enhance the framework. 

We presented here a quantitative analysis on the nature of noise for a differential equation model 
of prostate cancer, based on the construction of an agent-based version of the model. Specifically, we 
showed that the phenomenological interactions of this model are not able to explain observed noise in 
data. We suggested that a more detailed description of interaction and regulation mechanisms involved 
is needed to better clarify the noise effects. We plan to further investigate this direction, taking also into 
account spatial organization of the tumour. Our future work will also benefit from a comparison with 
experimental dataj^] 
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